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ABSTRACT 

In this paper we investigate whether Smoothed Particle Hydrodynamics (SPH), equipped with 
artificial conductivity, is able to capture the physics of density/energy discontinuities in the 
case of the so-called shearing layers test, a test for examining Kelvin-Helmholtz (KH) in- 
I— I stabilities. We can trace back each failure of SPH to show KH rolls to two causes: i) shock 

waves travelling in the simulation box and ii) particle clumping, or more generally, particle 
r K noise. The probable cause of shock waves is the Local Mixing Instability (LMI), previously 

. identified in the literature. Particle noise on the other hand is a problem because it introduces 

4^ a large error in the SPH momentum equation. 

Qh The shocks are hard to avoid in SPH simulations with initial density gradients because 

Q the most straightforward way of removing them, i.e. relaxing the initial conditions, is not vi- 

^ able. Indeed, by the time sufficient relaxing has taken place the density and energy gradients 

^ have become prohibitively wide. The particle disorder introduced by the relaxation is also a 

problem. We show that setting up initial conditions with a suitably smoothed density gradient 
dramatically improves results: shock waves are reduced whilst retaining relatively sharp gra- 
I dients and avoiding unnecessary particle disorder. Particle clumping is easy to overcome, the 

^ most straightforward method being the use of a suitable smoothing kernel with non-zero first 

central derivative. We present results to that effect using a new smoothing kernel: the Linear 
\Q Quartic (LIQ) kernel. 

We also investigate the role of artificial conductivity (AC). Although AC is necessary in 
^— I the simulations to avoid "oily" features in the gas due to artificial surface tension, we fail to 

\^ find any relation between using artificial conductivity and the appearance of seeded KH rolls. 

^ Including AC is necessary for the long-term behavior of the simulation (e.g. to get A = 1/2, 1 

KH rolls). In sensitive hydrodynamical simulations great care is however needed in selecting 
i-H the AC signal velocity, with the default formulation leading to too much energy diffusion. We 

^ present new signal velocities that lead to less diffusion. 

. The effects of the shock waves and of particle disorder become less important as the time- 

scale of the physical problem (for the shearing layers problem: lower density contrast and 
^ higher Mach numbers) decreases. At the resolution of current galaxy formation simulations 

mixing is probably not important. However, mixing could become crucial for next-generation 
simulations. 
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1 INTRODUCTION 

Smoothed Particle Hydrodynamics (SPH) is a lagrangian technique 
to solve the equations of hydrodynamics. It was first conceived by 
Lucy (1977) and Giiigold & Monaghan (1977) to solve astrophys- 
ical problems. Because the differential equations are solved on a 
particle mesh that moves with the flow, adaptivity to accomodate 
large density variations is inherent to SPH. As shown in Tasker 
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et al (2008) SPH indeed has the advantage over grid-based tech- 
niques when it comes to problems with a large dynamic range. 
However, when it comes to handling steep density gradients SPH 
codes have to acknowledge their superiors in grid-based codes. 

Recently, doubt has been cast upon whether SPH is able to 
capture all of the physics related to discontinuities. Agertz et al. 
(2007) highlighted a known problem in SPH: the formation of an 
artificial gap around density discontinuities. Previous allusions to 
this problem are found in e.g. ^ '" ^ ^ . • < (!^^^), where 
a method to enforce incompressibility in SPH is proposed. Tar- 
takovsky & Meakin (2005) propose a method to overcome this 
problem: use the local number density instead of the normal den- 



2 S. Valcke et ah 



sity to weigh SPH variables. This shifts the problem from requiring 
a continuous density to a continuous number density. The latter 
can be achieved by setting up the initial conditions such that the 
low density fluid uses SPH particles with lower mass, instead of 
less SPH particles with the same mass. Agertz et al. (2007) further- 
more argue that a basic SPH scheme is unable to exhibit Kelvin- 
Helmholtz or Rayleigh-Taylor instabilities when a density gradi- 
ent is involved. In a follow-up paper, Read et al. (2010) identify 
two main problems with standard SPH implementations: the Lo- 
cal Mixing Instability (LMI) and the "EO" error in the momen- 
tum equation. The LMI is the cause of the artificial gap problem 
highlighted above, whereas the problem with the EO error in the 
SPH momentum equation has previously been highlighted by 
ris (1996). Read et al. (2010) present a solution to these problems 
based on a temperature-weighted density and a modified smoothing 
kernel. 

Price (2008) presents a different solution to the artificial gap 
problem. Based on the general approach previously outlined by 
Monaghan (1997) he introduces an "artificial conductivity" (AC) 
term into the equations. This term induces a certain amount of en- 
ergy diffusion across energy discontinuities, allowing a disconti- 
nuity (which SPH is unable to handle properly) to become more 
"smeared out" and thus treatable with SPH. 

Kawata et al. (2009) present their implementation of an 
SPH scheme, closely tailored after the scheme by d & 
Price (2007), which includes the artificial conductivity term. They 
present the basic tests performed by Agertz et al. (2!i(/7): the blob 
test and the shearing layers test. Their results indicate that the 
Price (2008) AC solution gives improved performance for the blob 
test, on the condition that enough particles are present in the blob. 
Their shearing layers test (with a density contrast of 9.6) exhibits 
Kelvin-Helmoltz instabilities for a Kelvin-Helmholtz time-scale of 
TKH - 0.57. 

Interestingly, Okamoto et al. (2003) investigated shearing 
flows in SPH. They find that noisiness in the SPH smoothing of 
variables gives rise to small-scale pressure gradients which signifi- 
cantly decelerate the shearing flow. 

In § 2 we give a brief introduction to the SPH code, high- 
lighting the implementation of artificial conductivity (AC). A new 
smoothing kernel with non-zero central first derivative is presented 
in § 3. Various aspects of the shearing layers test are then examined 
in § 4 and a discussion is given in § 5. 

Most of the plots in this paper were made using HYPLOT. 
HYPLOT is a freely available^ open source analysis package, with 
an emphasis on SPH. Currently only GADGETII file reads are 
supported. HYPLOT uses PyQt4 for its GUI frontend, matplotlib 
for the plotting and a host of C++ classes for the actual computa- 
tions. Interested users are recommended to get the latest snapshot 
from the svn repository^. 



2 THE CODE 

We use the publicly available version of the Nbody/SPH code 
GADGETII (Springel 2005). There are two base premises when 
formulating a Smoothed Particle Hydrodynamics (SPH) solution 
to the equations of hydrodynamics: 

(i) the integral representation of field functions, 

^ http://sourceforge.net/apps/wordpress/hyplot/about/ 
^ https://hyplot. svn. sourceforge.net/svnroot/hyplot/trunk 



(ii) the particle approximation. 

In the first premise a function is replaced by its integral representa- 
tion, given by the integration of the multiplication of that function 
and a smoothing kernel function. The second premise states that 
the integral from the first premise is replaced by a discretized sum- 
mation using a set of particles in the support domain. The latter is a 
key approximation as it obsoletes the use of a background mesh for 
numerical integration. Both premises allow us to use the following 
simple equation to calculate the density at a certain point in space 
r: 

N 

p(r) = ^m,Vl^(|r-n|,/i), (1) 

i=i 

where the summation goes over all the particles within the support 
domain, delimited by the smoothing length h. Here, rrij is the mass 
of the jth particle, is a smoothing function. 

For a derivation of the SPH formulation of the basic equations 
of hydrodynamics we refer the reader to e.g. Springel & Hernquist 
(2002), where the equations are derived from a Lagrangian varia- 
tional principle. 

2.1 Artificial Conductivity 

From e.g. Monaghan (1997) and Price (2008) we learn that an arti- 
ficial conductivity (AC) term should be included in the SPH equa- 
tions when dealing with energy discontinuities. The expression for 
the dissipational part (i.e. without the adiabatic part) of the energy 
equation for an SPH particle i then becomes: 




+ OiuVsi^{ui - Uj)\^eij • ViWij, (2) 

with u the specific energy. The summation over j is the sum over all 
neighbours of particle i. Here, rrij is the neighbour particle mass, 
pij = {pi + pj)/2 with pi the standard SPH particle density, a and 
au are coefficients that can be used to dynamically vary the con- 
tribution of the terms based on the presence of e.g. local velocity 
convergence. We set a = = 1 and Vsig, = {ci-\-Cj — /3vij • Cij), 
with a the sound speed of particle i. Note that the first term in equa- 
tion (2) with 13 set to 1.5 is equal to the artificial viscosity (AV) 
employed in the GADGETII code (Springel 2005), apart from the 
Balsara switch. The second term in eq. (2) is the artificial conduc- 
tivity (AC) term. 

Price (2008) suggests the following form for the AC signal 
velocity Vs^^: 

^s-g = J^^, (3) 

V 

with Pi and Pj the pressures of respectively particles i and j. With 
this choice, spurious pressure gradients across contact discontinu- 
ities are gradually eliminated. We note however that as the expres- 
sion for the artificial conductivity ( ) is actually an SPH represen- 
tation of a diffusion term ( ), using the suggested signal 
velocity (eq. (3)) could lead to spurious energy diffusion. When 
using this signal velocity to eliminate pressure discontinuities one 
implicitly assumes that lower energy corresponds with lower pres- 
sure, as the sign of the energy transfer is determined by the energy 
gradient whilst the magnitude of the transfer is determined by both 
the energy and pressure gradients. If this assumption is valid the 
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diffusion will distribute energy from the high-energy particles to 
the low-energy particles, reducing the pressure discontinuities. In 
the simulations we use the equation of state for an ideal gas: 

P = (7 - (4) 

with p the pressure, 7 the adiabatic constant (typically taken to be 
5/3, the value for an ideal mono-atomic gas). From eq. (4) we learn 
that we can only know that higher energy corresponds to higher 
pressure when the density is constant. If the low-energy particles 
have higher pressure, the situation will not reach the desired pres- 
sure equilibrium as energy will flow from the high-energy particles 
to the low-energy particles, increasing the pressure discontinuities. 
In this case energy will keep flowing until an approximate energy 
equilibrium is reached. It is straightforward to formulate a signal 
velocity that does not suffer from this problem: 



:sign[(P.-P,)(t 



Pij 



(5) 



where the sign of the energy diffusion is now determined by the 
sign of the pressure gradient. Expanding on this idea, another pos- 
sible form of the signal velocity is: 



^sig,2 



: sign ^{Pi - Pj){ui - Uj)^ 



Pij 



where P is the SPH-averaged pressure: 



P. = ^mj^-pJ-^T^(|r-r.|,/i), 



(6) 



(7) 



with Aj the entropy of particle j, W = {Wi + Wj)/2. The ad- 
vantage of this form is that the artificial conductivity counteracts 
pressure differences at the fluid level, not at the particle level. A 
disadvantage is the extra overhead needed to compute and store P. 
Note that (as highlighted by Price (2008)) the formulations for the 
signal velocity shown here (eqs. (3), (5) and (6)) are not applicable 
when gravity is involved, because typically under hydrostatic equi- 
librium a configuration with constant pressure is not an equilibrium 
configuration. A signal velocity applicable in an environment with 
gravity, e.g. a modified version of eq. (3), can also be extended with 
the factor we added here in going from eq. (3) to eq. (5). 

It is worth considering whether the potential errors induced by 
equation (3) actually happen, and if so if their magnitude is of an 
order that requires fixing. To test this we set up two simulations, 
one using the signal speed of eq. (3), the other using eq. (5). These 
simulations consist of two central smoothed density gradients and 
constant pressure, exactly the same way we set up simulations for 
the shearing layers tests (section 4). Initial specific energies are set 
on a per-particle basis using the calculated SPH density and a fixed 
constant pressure value. To be able to clearly examine the influence 
of the conductivity the SPH particles are fixed in place. Results are 
shown in Fig. 1 . We immediately see that the simulation with the 
signal velocity as in eq. (3) (panel c) indeed exhibits divergent be- 
havior. SPH particles with low energy (0.25 < x < 0.75) receive 
extra energy, increasing their pressure, and vice versa for particles 
with high energy. Once this mechanism is set in motion it is self- 
amplifying. The simulation with the modified signal velocity on 
the other hand (panel d) does not exhibit any special behavior, the 
pressure remains constant, as expected. Note that the reason the di- 
vergent behavior is set in motion in this setup is numerical noise. 
Indeed, from eq. (3) one would expect a simulation with equal par- 
ticle pressures to have a signal velocity which is equal to zero ev- 
erywhere. As we set the initial particle specific energies based on 




Figure 1. a) Initial density along the y-axis for an x-slice [0,0.001]. b) 
Initial pressure along the y-axis. c) Pressure along the y-axis at time t = 2 
in the case of the standard AC signal velocity (eq. (3)). d) Pressure along 
the y-axis at time t = 2m the case of the modified AC signal velocity (eq. 
(5)). 



the required constant pressure value (P — 10) and the calculated 
density, we can expect that calculating back to the pressure later on 
(again using eq. (4)) will yield non-zero albeit very small pressure 
differences. This does not imply that the unwanted energy diffusion 
is a result of numerical artifacts, because in a real simulation there 
will always be non-zero pressure differences and we can expect at 
any time to have too much energy diffusion. On the other hand, be- 
cause in a real simulation the SPH particles are not fixed in place 
they will react to the increased pressure gradient, trying to erase 
it. The buildup of diffusion will thus not be as drastical as the one 
shown in Fig. I . We will investigate the new forms of the AC signal 
velocities further on. 



3 LIQ KERNEL 

A vital part of any SPH code is the smoothing kernel W. The at- 
tention given to the smoothing kernel is however somewhat limited, 
with people generally using the cubic spline kernel (see e.g. Kawata 
et al. 2009). Schuessler & Schmitt (1981) already demonstrated 
that the choice of the smoothing kernel can have a large impact 
on simulations. More recently, Morris (1996); Price (2005); Read 
et al. (2010) performed a linear perturbation analysis of the SPH 
equations of motion for different kernels, examining the stability of 
these kernels for longitudinal (related to the clumping instability) 
and transversal waves. They show several smoothing kernels with 
a zero central derivative that suffer from the clumping instability 
(also dubbed tensile instability). 

Various approaches have been suggested to deal with this. One 
possiblity is to modify the SPH equations. Monaghan (2000) advo- 
cates including an artificial pressure term. Sigalotti et al. (2009) on 
the other hand apply an adaptive kernel density estimation algo- 
rithm (ADKE). A different approach is to take a different form for 
the smoothing kernel. The advantage of the latter approach is that 
no modifications to the actual SPH scheme are necessary. The dis- 
advantage is that bell- shaped kernels (with central derivative equal 
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to zero) tend to give better results for a wide range of test cases as 
compared to hyperbolic or parabolic shaped kernels (Fulk 1996). 

Here we construct a new smoothing kernel with non-zero cen- 
tral derivative: the Linear Quartic or LIQ kernel. We note several 
earlier approaches: the HYDRA kernel by Thomas & Couchman 
(1992), who artificially modified the cubic spline (CS) kernel to 
have a constant central first derivative (by fixing it to — I/tt for 
X < 2/3). This approach was recently picked up by Merlin et al. 
(2009). Johnson et al. (1996) employed a quadratic kernel, result- 
ing in a linear first derivative. More recently Read et al. (2010) 
modified the central part of the cubic spline kernel with their Core- 
Triangle (CT) kernel, giving it a non-zero central first derivative. 
We will come back to these at the end of section 3.1. 

We assume a kernel of the form: 

Wr{u) 



W{r,h) = 



(8) 



with u — r /h and d the number of dimensions. For Wr we take the 
following functional form: 



< U < Xs 



Wr = Nx { f2 : Au^ + Bu^ + Cu^ + Du + E Xs < u < 1 







1 < u 



(9) 

Xs is a free parameter determining the connection point of the poly- 
nomial and linear functions. This form is inspired by two ideas: (i) 
we want the smoothing kernel to be smooth (ii) the first derivative 
of the smoothing kernel should be a monotonously ascending func- 
tion (i.e. it can have constant parts, but it can never descend). 

The equations for the smoothing kernel (9) have 7 free param- 
eters: A through F and the normalization factor N. We fix them by 
imposing the following boundary conditions: 



f2{Xs) 
df2. 



du 



df. 



—^{Xs) 



du 



Ml) 
"(1) 



(1) 



X 



W{r,h) dr 



du 






1, 



(10a) 
(10b) 

(10c) 
(lOd) 
(lOe) 

(lOf) 

(lOg) 



which ensure a smooth (up to second order) transition between /i 
and /2 at Xs and sufficiently smooth behavior for /2 ^ as ^ 1. 

Three different kernels are shown in Fig. 2 in the twodimen- 
sional case. The LIQ kernel has a lower central value than the CT 
kernel, and a higher value outwards. Read et al. (2010) find that 
they need a large number of neighbors for the CT kernel to reduce 
the Eq error of the momentum equation and to reduce the pressure 
blips at the boundaries. The Eq error on particle i is given by (Read 
et al. 2010): 



Pj 



hVW 



(11) 



Fig. 3 shows a plot of particle errors as a function of y, for the 
CS, CT and LIQ kernels. The number of neighbors is kept equal 
(32), the setup is that of the RH02 simulation (see Table 4 and 
Section 4). The difference between the CT and LIQ kernels on the 
one hand and the CS kernel on the other hand is big, with particle 
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Figure 2. Upper panel. The 2D LIQ, CS and CT smoothing kernels. Lower 
panel. First derivative of these kernels. 



clumping giving rise to large Eq errors, evidenced by the large scat- 
ter between the peaks. The results for the CT and LIQ kernels are 
comparable, with the LIQ kernel slightly in front, as can be seen 
from the reduced scatter away from the peak areas and the smaller 
maximum values overall. It would be interesting to perform the 3D 
simulations as performed by (Read et al. 2010) with the LIQ ker- 
nel, in this paper we restrict ourselves to 2D however. Both the LIQ 
and CT kernels have the desired properties for their respective first 
derivatives. Table 1 lists computed LIQ kernel coefficients for a 
range of Xs values. Analytical expressions for the coefficients can 
be found in appendix A. 



3.1 Fixing Xs 

To fix the remaining free parameter in the LIQ Kernel, the connec- 
tion point Xs, we set up a series of 2D Sod shock tube tests (Sod 
(1978), for recent use in SPH see e.g. Price (2008) (ID)). The x- 
range is [0 — 0.1], the y-range [0 — 1.5]. The adiabatic index 7 is 
set to 5/3. Further parameters can be found in table 2. Throughout 
this paper the number of SPH neighbours is 32, unless noted other- 
wise. Note that this number is actually quite high, with about 5-6 
neighbours per dimension. The particles are set up on rows in the 
^-direction i.e. the initial x-separation between particles is equal 
everywhere. Results are shown in figure 4. From these results we 
learn that varying Xs has a major impact on the simulation. A value 
of 0.3 is optimal, with higher values resulting in a more noisy den- 
sity profile and smaller values leading to small density oscillations 
in certain parts of the profile. The shock tube result with the LIQ 
kernel and Xs = 0.3 is slightly less good than that for the CS kernel 
(bottom right panel). The density found by the SPH summation us- 
ing the LIQ kernel is less accurate than the density found using the 
CS kernel. The simulations in figure 4 started from the same initial 
conditions file, but have, depending on the value of Xg, slightly dif- 
ferent starting densities: the density is overestimated, the more so 
for smaller values of Xg. Instead of lowering the particle masses to 
get identical initial density profiles we use identical masses in all 
simulations. The analytical shock-tube results shown in Fig. 4 are 
computed using the actual initial SPH densities. 

In Fig. 5 we show the particle distribution on top of the ren- 
dered density, for a small region of the shock tube, at t = 0.2. 
Overplotted is the circle of interaction for the same SPH particle. 
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Table 1. Coefficients of the LIQ Kernel (eq. (9)). The last two columns (A^2D, A^3d) 
give the norm respectively for two and three dimensions. See appendix A for 
analytical expressions. 



Xs 


A 


B 


C 


D 


E 
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A^2D 


A^3D 
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1 





-1 


05 


0.5 


4.775 
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0.9259 


0.4630 


0.7 
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0.8 
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Figure 3. The Eq error for three shearing layers simulations (t = 0.5) with, 
from top to bottom: the cubic spline kernel, the core triangle kernel, the 
linear quartic kernel. To guide the eye a binned mean value is also shown. 
The cyan line, shown in all panels, shows the values for the CS kernel. The 
green Hnes in panels b and c show the values for the respective CT and LIQ 
kernels. 



Table 2. Sod shock tube test parameters. (1) density (2) 
specific energy (3) pressure (4) number of particles per y- 
column (5) total number of particles 
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P 


A^col 


N 


high-density 
low density 
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Figure 4. Results for the 2D Sod shock tube test with the LIQ Kernel, using 
varying connection points (xs). Only particles in the x-interval [0,0.005] are 
shown (x-range: [0,0.1]). Connection point values are shown in the plots. 
Panel h) shows the result using the standard Cubic Spline kernel. The over- 
plotted solid green line in each plot is the analytical solution. 



The clumping behavior in the left plot, using the CS kernel, is strik- 
ing: particles form groups of 2. This behavior has been known for 
some time (see e.g. Schuessler & Schmitt 1981). When using the 
LIQ kernel (right panel) there is no clumping. We note that this 
clumping does not lead to an actual decrease in resolution, as the 
circle of interaction for the SPH particle has a comparable radius 
in the two plots. Indeed, the clumping does not affect the actual 
size of the smoothing region. It does affect the homogeneity of the 
particles in that region, and through that it affects the validity of 
using an SPH estimate of a continuous integral (see section 2). As 
the CS kernel has been used extensively in SPH codes and found to 
give excellent results (see e.g. Springel 2005; Price 2008; Kawata 



et al. 2009, the shock tube results in Fig. 4), the SPH estimate of 
the integral does not seem to suffer from the clumping. 

In Fig. 6 we show the energy conservation of SPH for the 2D 
Sod shock-tube test for the HYDRA, Core Triangle, Johnson, LIQ 
(xs = 0.3) and Cubic Spline kernels. As expected, the Cubic Spline 
kernel clearly leads to the best energy conservation. The behavior 
of the CT and LIQ kernels is similar, especially in their initial evo- 
lution. The Johnson and HYDRA kernels show the largest growth 
in their respective relative energy errors. Note that the HYDRA ker- 
nel is identical to the CS kernel, apart from the modified central first 
derivative. The resulting kernel does prevent clumping but the en- 
ergy conservation leaves much to be desired. Moreover, among the 
family of LIQ kernels, the choice Xs = 0.3 turned out to provide 
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Figure 5. Particle distribution on top of the rendered density for a region 
in the shock tube test. Left panel: Cubic SpHne. Right panel: LIQ Kernel. 
The overplotted circle is the region of influence for the same particle in both 
runs. The particle clumping in the left panel, using the CS kernel, is striking. 
No clumping is seen in the right panel. 
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Figure 6. Relative energy error as a function of time, for a number of 2D 
shock-tube tests, using different smoothing kernels. Lines indicate, from top 
to bottom: HYDRA kernel (CS with modified central derivative), Johnson 
kernel (quadratic), LIQ kernel (see section 3, linear + 4* order), CT kernel 
(centrally modified CS), Cubic Spline kernel (3^^ order). Note how modi- 
fying the central value of the derivative of the CS kernel drastically reduces 
energy conservation. 



the best energy conservation. The different result for the LIQ and 
CT kernels arises between t = 0.07 and t = 0.15. This interval 
corresponds to the time where the initial rectangular symmetry in 
the simulation breaks up. The initial energy error (t < 0.07) thus 
shows the performance of the SPH scheme on a regular particle dis- 
tribution whilst the late energy error (t > 0.15) shows the perfor- 
mance with SPH particles having arranged themselves according 
to the respective smoothing kernels. The poor performance of the 
HYDRA kernel can be attributed to the artificial modification of the 
central derivative of the kernel, which breaks the conservative form 
of SPH. Note that all simulations use individual time-steps which 
in itself breaks the time- symmetry (and hence the energy conserva- 
tion) of the integrator. 



4 SHEARING LAYERS TEST 

The shearing layers test, as presented by Agertz et al. (2007), is an 
excellent way of examining the ability of a hydrodynamical code 



to capture the formation of Kelvin-Helmholtz instabilities. These 
instabilities occur when two fluids have a relative velocity perpen- 
dicular to their contact layer. For an incompressible fluid there is 
an analytical expression for the time-scale for the growth of these 
Kelvin-Helmholtz instabilities (Chandrasekhar 1961; Price 2008): 

where 



tkh 



(12) 



27r (pi/>2)2 \Vx,l - Vx,2\ 



A 



{pi +P2) 



(13) 



The indices 1 and 2 denote the two fluid layers. A is the wavelength 
of the perturbation. When re- written in terms of the density contrast 
X (pi — XP2, we make the convention x ^ 1 ^-g- the index 1 
denotes the high-density layer): 



271 (x)2 \Vx,l - 



A 



(1 + x) 



(14) 



Various approaches have been suggested to enforce incompressibil- 
ity in SPH (e.g. Cummins & Rudman 1999; Hu & Adams 2009). 
We use the weakly compressible SPH formulation, approximating 
incompressibility by using a small Mach number. Note that incom- 
pressibility is a requirement only for equation (14) to be valid, it is 
not a requirement for the simulations themselves to be valid. 

KH instabilities are seeded by introducing an initial ^/-velocity 
perturbation following the prescription: 



Asin(27rx/A). 



(15) 



All simulations, unless noted otherwise, have the following setup: 
the X- and y ranges are 1 . For the adiabatic index 7 we use a value 
of 5/3, as applicable for an ideal monoatomic gas. The density con- 
trast X is 10, with low- and high-density layers having a density 
of respectively 1 and 10. The respective specific energies are 15 
and 1.5, resulting in respective sound speeds of 4.08 and 1.29. A is 
fixed at 1/6, which should result in the initial growth of 6 KH den- 
sity rolls. A is set to 0.025. Note that we focus on a density contrast 
of 10, not 2, because preliminary tests showed us that an SPH code 
has much less trouble reproducing Kelvin-Helmholtz instabilities 
for lower density contrast: we concentrate on the more demand- 
ing test here. Previous studies of the KH instabilities either restrict 
themselves to a density contrast of 2 (Read et al. 2010; Merlin et al. 
2009) or show KH tests with KH time- scales of 0.6 or lower (Price 
2008;Kawata et al. 2009). 



4.1 Grid Code 

To get a reference to compare the SPH results with we have con- 
ducted a series of Eulerian hydrodynamical simulations, using the 
FLASH code (version 3.2, (Fryxell et al. (2000); Dubey et al. 
(2008), see also http://flash.uchicago.edu)) . FLASH 
is an adaptive mesh refinement hydrodynamics code. Its standard 
hydro- solver uses the piecewise-parabolic method. 

The KH simulation is run on a periodic 2D grid of size 1x1. 
In order to avoid confusion from numerical issues as much as pos- 
sible, we switched off the steepening algorithm for contact discon- 
tinuities and use just one refinement level, resulting in an effective 
resolution of 256^ = 65536 grid cells. We note that differences to 
runs with 2 refinement levels and contact steepening appear only 
after a few KH timescales, but are evident in the late stages of the 
higher Mach number runs presented here. The initial conditions are 
set as in the SPH simulations, without taking special care to smooth 
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Table 3. Velocities and tkh for the different KH sim- 
ulations. Rows from top to bottom: suffix: simulation 
number suffix (e.g. for the SPH series: SPH1-SPH5), 
M: Mach number of the high-density layer, \vx\'. ab- 
solute value of the x-velocity, tkh- KH time-scale 
(eq. (12)) 



suffix 


1 


2 


3 


4 


5 


M 


0.2 


0.4 


0.6 


0.8 


1 




0.26 


0.52 


0.77 


1.0 


1.3 




1.23 


0.56 


0.37 


0.28 


0.22 



gradients, but keeping the idealised sharp discontinuity in density 
and velocity. Results are shown in Figs. 6 and 7 for two simulations 
with respective Mach numbers of 0.2 and 0.6. Those KH rolls that 
appear do so at times consistent with their respective tkh. 

Results are shown in Figs. 7 and 8 for two simulations with re- 
spective Mach numbers of 0.2 and 0.6. Those KH rolls that appear 
do so at times consistent with their respective tkh . 



4.2 Standard SPH 

We first examine the capability of standard SPH to handle the shear- 
ing layers problem. With standard SPH we mean SPH as imple- 
mented by e.g. Springel (2005) (which uses the entropy formula- 
tion of SPH) augmented with the artificial conductivity formalism 
(see section 2). This form of the SPH equations is rapidly becoming 
the standard (see e.g. Price 2008; Kawata et al. 2009; Vanaverbeke 
et al. 2009; Merlin et al. 2009). We do this test for a range of Mach 
numbers, and thus for a range of Kelvin-Helmholtz time-scales (eq. 
(12)). The specific energies of the SPH particles were set after the 
density is computed, in order to force perfect pressure equilibrium 
(thus applying a crude smoothing of the IC). Details for the differ- 
ent simulations are listed in tables 3 and 4. Rendered plots of the 
density are shown in Fig. 10, for each run at its respective tkh (see 
table 4). From Fig. 10 we learn that changing the Mach number in 
the simulation has a drastic impact on the formation of KH rolls: 
for small M they are completely absent. The simulations of Price 
(2008) with X = 10 have M = 0.775 for the high-density layer, 
they are to be compared with SPH4. We thus find that the results 
in Price (2008), where artificial conductivity was enough to enable 
SPH to produce KH rolls, are only valid for high enough Mach 
numbers. In following sections we will try to improve on this situa- 
tion. The situation at t = 2 is shown in Fig. 1 1. It tells us that there 
is a large difference in the long-term evolution of the shearing layer 
simulations depending on the setup and the ingredients of the code. 
On the panels for M = 0.6 (c, h and m) we should see A = 1/2 
instabilities when comparing with the grid simulations in Fig. 8. 
Panel h shows a hint of these instabilities but they are very poorly 
developed. Panel m does better, the KH rolls appear although only 
on one side of the high-density layer. 



4.3 Standard SPH - Column-smoothed IC 

We can expect the way the initial conditions for the shearing layers 
test were set up in the previous section to cause problems (as briefly 
mentioned by Hess & Springel (^^^^)^ though they do not try to for- 
mulate a remedy). SPH is inherently smoothed, and, in the current 
formulation, this particle smoothing is isotropic as the smoothing 
length is a scalar, not a vector nor a tensor. When one uses this 



Table 4. Parameters for the different KH simulations, series: Name for 
the simulation series, e.g. SPH1-SPH5. #part'. Number of particles used. 
AC: Artificial conductivity included? (yes: X). kernel: Smoothing kernel 
used, p-smoothing: Type of density smoothing applied, if any. 



series 


SPH 


RHO 


LIQ 


GRID 


NOAC 


#part 


199252 


184100 


184100 


200788 


184100 


AC 


X 


X 


X 


X 




kernel 


CS 


CS 


LIQ 


CS/LIQ 


LIQ 


p-smoothing 




column 


column 


grid 


column 



formulation of SPH to tackle a problem where sharp discontinu- 
ities are present, as is the case in the shearing layers test, one is 
forcing SPH into a situation where its behavior is not well defined. 
Smoothing the initial particle energy based on the calculated initial 
density and a requested constant pressure value does not remedy 
this problem: it enforces equal pressures on a particle basis, which 
does not necessarily result in pressure equilibrium on a simulation 
basis. This is easy to see when making an analogy with the density: 
giving all particles an equal mass results in a constant density only 
if the particle distribution is homogeneous. Robertson et al. (2010) 
show that Eulerian codes exhibit similar problems, with numeri- 
cal diffusion wiping out small-scale structures in the presence of a 
large bulk flow and sharp discontinuities. They present a different 
version of the KH test which does not suffer from this. We however 
do not discard the canonical test because the code is not able to 
solve it. We choose to modify the test so that we retain the physics 
of the original test whilst allowing the code to deal with it 

Artificial conductivity goes a long way in preventing discon- 
tinuities as sharp as these from arising in simulations, but by the 
time AC has sufficiently smoothed the initial discontinuity a se- 
ries of shocks has been created in the simulation volume. These are 
shown in Fig. 1 Depending on the time-scale of the growth of the 
instabilities these shocks will play a role in destroying emerging 
KH instabilities: the longer it takes for the instabilities to manifest 
themselves, the more they are wiped out by the travelling shock 
waves. When the shocks pass through the emerging instabilities a 
mixed layer is formed at the contact discontinuity. This layer acts 
as a lubricant, separating the two layers and thus preventing KH 
instabilities from developing. Only in the case where the KH in- 
stabilities form fast enough compared to the scale on which they 
are destroyed are they able to survive. This can be translated into 
a general paradigm: the longer it takes for KH rolls to manifest 
themselves in SPH, the more time there has been for SPH-induced 
deviations from the analytical, theoretical problem to destroy them. 
This explains why it is more easy to produce KH rolls for high 
Mach numbers. The origin of these shocks is probably related to 
the LMI, in the same way as the artificial gap problem. Imagine two 
particles with equal pressures, but different densities. If these par- 
ticles approach each-other their respective densities will change to- 
wards equality. As entropy is conserved this implies, with^ = Ap^ 
(Springel & Hernquist 2002), that these particles now have differ- 
ent pressures. As this process happens everywhere along the initial 
contact discontinuity this can set a shock-wave in motion. A pos- 
sible cure for this problem is to smooth the regions where these 
discontinuities occur, giving rise to a net smaller LMI effect. To 
create a smooth interface we choose a function of the form: 

p{y) = A atan(^(y + C)) + D, (16) 
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Figure 7. Grid simulation of the shearing layers test with M = 0.2. Top row, left to right: t = 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2. Bottom row, left to right 
t = 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4: White boxes mark blocks of 16^ grid cells. 
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Figure 8. Grid simulation of the shearing layers test with M = 0.6. Top row, left to right: t = 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2. Bottom row, left to right 
t = 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4: White boxes mark blocks of 16^ grid cells. 



where 



A = 


pi - po 


(17) 




2 atan(/3) 


B = 




(18) 


C = 


,5 
2 


(19) 


D = 


pi + P2 

2 


(20) 



were fixed requiring (16) to be symmetrical around the boundary, 
scaling the argument of the at an such that y = 0,6 resulted in 
an argument of respectively —/3, 13. For /3 a value of 10 gives rea- 
sonable results, because atan(lO) is close to 7r/2 (which is needed 
because it is connected to two constant-density functions). 6 is the 
width of the boundary, which we took to be (^/max — Vmin) /7.5 for 
a part of the simulation box containing a single boundary layer e.g. 
S = 0.5/7.5 in the current setup. 

There are several methods to set up a smooth density. One 
method is to set up the particles in the two regions on two differ- 
ent grids. Within such a grid the x- and y separations (Ax and 
respectively) of particles are equal. Denoting one layer with the 
subscript 1, the other with 2, we have of course Axi / Ax2 and 
A^i / Ay2 when the layer densities differ. The smooth boundary 
between those grids can then be constructed by putting a number of 
equal-^ rows between those grids. The distance between these rows 
varies from Ayi to Ay2 according to the square (in the 2D case) of 
some chosen function (e.g. eq. (16)). The number of particles on a 
row is then fixed by comparing the computed SPH density of par- 
ticles on that row to the value of required analytical density pan at 



that point. To get a perfect representation of the bounding function 
one would have to use an iterative scheme, varying the number of 
particles on the rows, recomputing the smoothing lengths and den- 
sities and repeating until convergence is reached. In our setup our 
prime concern is that the density is smooth, it is of little importance 
if the analytical boundary function is perfectly reproduced. We thus 
use a shortcut in generating the initial conditions: the SPH density 
on the i-ih row pi is esimated by: 

(21) 



AxiAyi 



with rrii the particle mass on that row, Ayi the distance to the pre- 
vious row and Axi the interparticle distance. Axi, and from that 
the number of particles, can then be found from: 



Axi 



pan Ayi 



(22) 



Another method to construct a smooth density is to put all the 
SPH particles on columns (lines with constant x), with an equal 
and unchanged intercolumn distance. The smoothed density can 
then be attained by varying the interparticle separation within those 
columns. Finding those separations, using equation 21, is then a 
simple matter of using some bisection algorithm. 

Simulations RHOl-5 correspond to the SPHl-5 series, with 
their initial density discontinuity smoothed using the first smoothed 
density method: particles on columns. Results are shown in Fig. 10. 
When comparing the top and middle rows we see that for M > 0.6 
KH rolls have appeared. For M = 0.4 some bumps are present, but 
these do not evolve into KH rolls. Although we have already greatly 
improved the situation, for low Mach numbers we are still unable 
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Initial specific energy (a), density (b) and pressure (c) as a function of y for the SPHl-5 (see table 3) simulations. Only particles in the x-interval 
are shown. Note that the specific energy of the particles at the contact discontinuity is smoothed in order to avoid a pressure blip. 
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Figure 10. Density plots for shearing layers simulations (table 4) at their respective tkh (table 3). Top row (a-e): from left to right: simulations SPH1-SPH5: 
standard SPH simulations. Middle row (f-j): RH01-RH05: standard SPH with a smoothed density setup. Bottom row (k-o): LIQl-5: SPH with a smoothed 
density and using the LIQ kernel. 



to get adequate results. In Fig. 12 we can see that the magnitude of 
the shocks has decreased because of the density smoothing. 

4.4 Modified kernel - Column-smoothed IC 

As mentioned in section 3 we can expect to do better in very sensi- 
tive simulations when using a smoothing kernel that does not cause 
particle clumping. In Fig. 10 we show results for the same simu- 
lations as the RH01-RH05 simulations (e.g. with smoothed den- 
sity), but now using the LIQ kernel (see section 3). We see KH 
instabilities appear for M > 0.4, which is an improvement over 



the previous results where clear KH rolls appeared for M > 0.6. 
We also see that in general the shape of the KH rolls is much more 
symmetrical. 

In figure 13 we show the pressure of the particles at roughly 
half the KH time- scale. When we look at the column smoothing 
results it is clear that the dominant shock in the RH02 simulation, 
around y = 0.5, is greatly reduced in magnitude in the LIQ2 simu- 
lation. This difference is mainly due to the different way the kernels 
respond to the initial conditions (particles initially on columns): 
the CS kernel is clearly not a good choice in this case. This is 
also demonstrated in the bottom right panel of Fig. 13, where a re- 
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Figure 11. Same figure as Fig. 10, all figures now att = 2. 
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Figure 12. panel a. simulation SPHl at t = 0.1 (particles limited to x-interval [0,0.005]). One clearly sees shocks travelling inward from the contact 
discontinuities, panel b. simulation RHOl. Shocks are still present, although reduced in magnitude, panel c. density smoothing + the LIQ kernel. 

suit for the CS kernel using grid smoothing is shown. It compares defined. The situation for LIQ+1 however has not changed: KH 
favourably with both LIQ2 results. rolls are still abscent at that resolution. 



4.4.1 Particle number 

We can try increasing the particle number for the simulations with 
low Mach number, to see how it affects the KH rolls: simulations 
LIQ+1 and LIQ+2, each with 551 100 particles. Results are shown 
in Fig. 14. The LIQ+2 simulation, where KH rolls were already 
present in its low-resolution counterpart LIQ2, is seen to benefit 
a lot from the increased resolution: the KH rolls are much better 



4.5 Grid-based density smoothing 

We have also constructed initial conditions using the second 
smooth interface method: the two different layers consist of two 
different grids, with two smooth interfaces at the contacts. The en- 
ergies of the particles are set using the same algorithm as before: 
after the SPH densities are computed the particle specific energies 
are set based on a constant pressure value. Part of the initial par- 
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Figure 15. Initial conditions with smooth boundary from section 4.5 Left 
panel: At the top and bottom we can see the high- and low density grids 
respectively. A smoothly varying layer lies in between. Right panel: SPH- 
computed density as a function of y, for the particles in the left panel. The 
density varies smoothly. 
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Figure 13. Pressure of SPH particles in the x-interval [0,0.005] as a func- 
tion of y, diit = 0.25, with M = 0.4. The top row shows two runs with 
LIQ2 setup. The bottom row shows RH02. The left panels have column 
density smoothing, the right panels have grid-based density smoothing. 
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Figure 14. Density of two increased resolution shearing layers simulations 
at their respective tkh- panel a. LIQ+1 (M = 0.2, 551 100 particles). 
panel b. LlQ+2 (M = 0.4, 551 100 particles). 



tide distribution is shown in Fig. 15. The simulation shown uses 
the LIQ kernel. Results for M = 0.2 are shown in Fig. 16. On the 
density plot we see that no KH rolls manifest themselves: the same 
shocks that were present in the LIQl simulation are present here. 
We thus find that the artificial conductivity is not able to smooth 
away initial energy discontinuities fast enough to prevent the LMI 
from triggering these shock-waves. In the right panel we show the 
rendered density of the RH02 simulation, using grid-based density 
smoothing (see also Fig. 13). The observed ripples at the contact 
discontinuities are a slight improvement over the RH02 simulation 
shown in figure 10. The particle clumping does however prevent the 
KH rolls from growing to the same size as their LIQ2 counterparts. 

In figure 1 7 we show the destruction of the KH instabilities by 
a shock wave in detail. From bottom to top we can see the wave 
passing through the emerging KH instabilities, erasing them. 




X X 



Figure 16. Shearing layer simulations at their respective tkh using grid- 
based density smoothing. Left panel.M = 0.2, using the LIQ kernel. Right 
panel: M = 0.4, using the CS kernel. No improvement is seen over the 
column-based smoothing results. 



4.5.1 Relaxing 

Because we want to avoid any influence from shock waves, which 
we have shown are moving through the simulation box, we set up 
a simulation where we have first relaxed the initial conditions. The 
relaxing scheme we used is a default SPH simulation, with the PdV 
terms removed from the energy equation. The energy of particles is 
thus not allowed to change due to expansion/contraction, only due 
to artificial viscosity and conductivity. We use this scheme to stay 
as close as possible to the original energy profile. 

We have relaxed run LIQl. Two different runs were started 
from this relaxing run, the first starting from the relaxing run snap- 
shot at t = 0.5, the second from the t = 1.0 snapshot. None of 
these simulations show any formation of KH rolls. Rendered views 
of Vy are shown in respectively Figs. 18 and 19. In Fig 18. we can 
see that the magnitude of the shock (bottom panel, around y = 0.7) 
is greatly reduced compared to the shock in Fig. 17. The emerging 
instabilities are also reduced in magnitude when compared to those 
in Fig. 17, due to the increased particle disorder and widening of 
the discontinuities. The small shock is sufficient to remove these 
smaller emerging instabilities. In Fig. 19 we can see that the mag- 
nitude of the KH seeds is very small. They do not develop and 
slowly fade away. This is the result of the long relaxing time which 
greatly diffuses the energy discontinuity and introduces too much 
particle disorder. 
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Figure 20. Shearing layers simulations without artificial conductivity (column-smoothed density). Left to right: Simulations NOAC1-NOAC5. Top row: 
simulations at their respective tkh- Bottom row: simulations at t = 2. 



4.6 Artificial Conductivity? 

We have tested the need for artificial conductivity to be included in 
shearing layers simulations using the LIQ kernel and with an ini- 
tially smoothed setup. To do this we repeated the LIQ 1-5 simula- 
tions, with AC switched off. Results are shown in Fig. 20. The com- 
parison with the bottom row of figure 10 learns us several things: 
i) applying artificial conductivity does not lead to the appearance 
of KH rolls. It is possible that with initial conditions which are not 
smooth and when using the CS kernel AC, provides an amount of 
smoothing in time for KHs to develop where they would not oth- 
erwise, ii) Artificial conductivity is needed to avoid an "oily" na- 
ture of the gas. iii) When comparing the long-term evolution of the 
simulations (Figs. 1 1 and 20), artificial conductivity in its default 
implementation is both a blessing and a curse. The panels in the 
bottom row of Fig. 20 show small-scale structures: "holes", inclu- 
sions of low-density gas inside the high-density layer, are present 
in the high-density medium. When comparing with the grid-based 
simulations (Figs. 7 and 8) we see that these are also present there. 
In the simulation with AC applied these holes have entirely disap- 
peared. It is however obvious that the tendency for the layers to 
avoid mixing without AC prevents realistic long-term behavior. 



4.6.1 Gap? 

As recently highlighted by Agertz et al. (2007), SPH in its stan- 
dard form (without AC) suffers from the formation of a "gap", i.e. 
a small void layer between two layers with different densities. In 
Fig. 21 we show plots of two simulations without artificial conduc- 
tivity, one with the CS kernel, the other with the LIQ kernel. It is 
clear that the LIQ kernel prevents the formation of a wide gap. As 
previously shown there still is a need for artificial conductivity (see 
Fig. 20) because the layers are still reluctant to mix. As we already 
achieved a great deal of improvement by trying one new smoothing 
kernel, we do not think it impossible that still better formulations 
are possible. 



4.6.2 Signal velocity 



In section 2. 1 we briefly discussed the signal velocity used when in- 
cluding AC in simulations. In Fig. 22 we show some results using 
the modified signal velocities for simulations with M = 0.6 and 
LIQ3 setup. They should thus be compared to panel m of Fig. 11. 
The top two rows show simulations using ^^^ig,i (eq. 5), the bottom 
two rows show results using ^;^ig,2 (^Q- 6). For both signal veloci- 
ties, the results for t = tkh are in good agreement with the grid 
simulation results as well as the other SPH results. At t = 2 the 
results differ from the grid results, with the A = 1/2 clearly sur- 
facing in the SPH simulations whereas this instability takes longer 
to surface in the grid results. When comparing with Fig. 11 it is ob- 
vious that for both signal velocities the amount of energy diffusion 
applied is smaller. Indeed, there is much more small-scale structure 
left in the t = 2 panels in Fig. 22. We show a series of snapshots at 
different times in Fig. 23, in order to have a clear view of which in- 
stabilities are surfacing when. The top two rows (panels a-p) show 
the evolution of the LIQ3 simulation. It goes from A = 1/6 to 
A = 1/2 rolls, although the latter only appear on one side of the 
high-density layer. Panels q-F show a LIQ3 simulation with ^^^ig,i 
AC signal velocity (eq. (5)). Changing the sign of the signal ve- 
locity is a clear improvement in this case, with the KH rolls being 
even better resolved than in the LIQ3 simulation and with the rolls 
appearing on both sides of the central flow (see the grid simulation 
results in Fig. 8). Panels G-V show a LIQ3 simulation with ^^^ig,2 
signal velocity (eq. (6)). The KH rolls at t = 2 are less well defined 
than in the LIQ3 simulation, they do appear on both sides of the 
flow. The result using ^^^ig,i thus compares favorably to the result 
using i^^ig,2 and using v^^g^i requires virtually no extra computa- 
tions, as opposed to using Vsig^2- It seems that using v^ig^i allows 
AC to act sufficiently strong to prevent clear layer separation ("oili- 
ness" due to the LMI), but not too strong so as not to lose too much 
resolution by diffusing too much energy. Note that Fig. 22 shows 
that the actual results can change drastically when using a different 
resolution. 
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Figure 23. 16 snapshots of the density of 3 shearing layers simulations. Top rows (a-h, q-x, G-N): t = 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, bottom rows (i-p, 
y-F, O-V): t = 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4. panels a-p: simulation LIQ3. panels q-F: simulation LIQ3 with modified AC signal velocity v^-^^ (see 
also Fig. 22, vilowres). This simulation develops clearly resolved A = 1/2 KH rolls, panels G-V: simulation LIQ3 with modified AC signal velocity v'^^^ ^ 
(see also Fig. 22, V2lowres). This develops less well defined A = 1/2 KH rolls, these do appear on both sides of the central stream, contrary to the LIQ3 
simulation. 



5 DISCUSSION 

Using a large suite of simulations we have shown that the perfor- 
mance of SPH on the shearing layers test can be improved drasti- 
cally, if i) the density gradient in the initial conditions is smoothed 
and ii) a smoothing kernel that does not cause particle clumping is 
used. Two effects are in play that can cause things to go haywire: 
shocks travelling through the simulation box and particle clump- 
ing, or more general, particle disorder. The effect of SPH particle 
disorder was already discovered by Okamoto et al. (2003), who 
set up a simulation with a small, hot shearing flow layer inside a 
cold medium. They found that noisiness in the SPH smoothing of 
variables gives rise to small-scale pressure gradients which signif- 
icantly decelerated the shearing flow. As using the Cubic Spline 
smoothing kernel causes particles to clump together in groups of 2, 
as can be seen in the left panel of Fig. 21, we can expect deceler- 
ation of SPH particles located in a small layer around the contact 
interface. This layer then acts as a lubricant between the two shear- 
ing layers, removing direct contact of the latter two and hence pre- 
venting KH instabilities from forming or growing. The LIQ kernel 
on the other hand gives rise to a much more homogeneous particle 



distribution (right panel of Fig. 21). Employing a suitble smooth- 
ing kernel with non-zero central first derivative (we say suitable 
because although the first derivative being non-zero is necessary 
to avoid clumping, it is not sufficient: a LIQ kernel with connec- 
tion point Xs = 0.5 suffers from even heavier clumping than the 
CS kernel: particles tend to clump together in groups of 6) gives 
rise to a much more homogeneous particle distribution. This has 
a major impact on shearing layers simulations, allowing KHIs to 
form much easier in general. The effect of particle disorder has 
been highlighted in the literature on various occasions. Read et al. 
(2010) show that particle clumping gives a large EO error in the 
SPH momentum equation, which in turn prevents mixing in SPH. 
Morris (1996) also found that that SPH does not converge at flow 
boundaries because of the EO error. 

The shock problem is not so easily tackled. The shocks are 
triggered by the local mixing instability, which occurs because en- 
ergy (entropy) is not smoothed in SPH. Several authors have iden- 
tified this problem (e.g. Cummins & Rudman 1999; Tartakovsky & 
Meakin 2005; Agertz et al. 2007; Price 2008; Read et al. 2010) and 
various solutions have been suggested. We examined the results us- 
ing the artificial conductivity solution (Price 2008). The magnitude 
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Figure 17. Rendered Vy for a simulation with M = 0.2, LIQ kernel, grid- 
based density smoothing. From bottom to top: t = 0.35, 0.4, 0.45. Bot- 
tom panel: We see emerging KH instabilities as the 6 red Vy areas around 
y = 0.75. Below them, around y = 0.7, we have a shock wave travel- 
ling upwards, with a small layer of positive Vy in front of it and a region 
of negative Vy in its wake. Middle panel: The shock is now almost at the 
same height of the instabilities. Its front is entering the low-density layer, 
increasing its y-length. Upper panel: The shock has passed through the in- 
stabilities. These now have a drastically reduced magnitude. The wake of 
the shock is still visible as the big blue band at the top of the panel. 



of the shock-waves can be reduced by applying an initial density 
smoothing, which reduces the magnitude of any shock waves that 
might develop, leading to a significant increase in the simulation 
results. To get rid of the remaining shocks one can relax a simula- 
tion before applying velocities to it. This does not help in this case 
however because relaxing not only significantly widens the density 
and energy discontinuities, it also removes a great deal of the sym- 
metry that was initially present in the problem, i.e. it introduces a 
lot of particle disorder. Indeed, a simulation started from a grid has 
perfect initial particle symmetry and as sharp a discontinuity as one 
desires. When relaxing the initial conditions on the other hand, par- 
ticle positions are shifted to reach an equilibrium. The smoothing 
kernel plays an important part in this because through small-scale 
variations in forces it will determine the final configuration of the 
particles. This can readily be seen in Fig. 5. Even when using the 
LIQ kernel there is an increased amount of particle noise. When 
bulk shearing layers velocities are then given to particles based on 
on which side of a line they lay, and velocity perturbations applied 
to particles irrespective of the underlying configuration of the par- 
ticles, it is not hard to imagine that the problem of the particle noise 
will quickly result in momentum transfer at the contact layers, thus 
shutting down all KH-related activities. Furthermore, we found that 
the onset of the shock-waves, occuring because of the LMI trig- 
gered in the initial particle configuration, is not removed by the 
artificial conductivity. 

We found that adding artificial conductivity to the SPH 
scheme does not have an impact on whether or not initial ( A = 1/6) 
KHIs surface. Although it has been reported and used as such in 



Figure 18. M = 0.2 simulation starting from a LIQl IC which was relaxed 
during a time of 0.5. From bottom to top: t = 0.35, 0.4, 0.45 (after the end 
of relaxation). The same shock wave seen in Fig. 17 is present. However, its 
magnitude is greatly reduced, due to the relaxing of the IC. The initial KH 
seeds are also reduced in magnitude due to the increased particle disorder 
and general widening of the discontinuities. 
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Figure 19. M = 0.2 simulation starting from a LIQl IC which was relaxed 
during a time of 1. Bottom panel: t = 0.35 (after the end of relaxation). 
Top panel: t = 0.5 (after relaxation). No shock is seen, the initial particle 
disorder and widening of the discontinuities themselves remove the growing 
KH instabilities. 



the literature (see e.g. Price 2008; Kawata et al. 2009) we show that 
simulations that use a suitable smoothing kernel and a smoothed 
initial density gradient are equally able to form these KHIs, ir- 
respective of AC being included or not. Even the formation of a 
visible "gap" (Agertz et al. 2007) is prevented by using a suitable 
smoothing kernel (Fig. 21). Including AC is, with the smoothing 
kernels used in this paper, still necessary to i) allow mixing to hap- 
pen, avoiding "oily" features in the SPH gas phases (see Fig. 20) 
and preventing the LMI from triggering during the course of the 
simulation and ii) get the A = 1/2 KH rolls later on. Here, point 
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Figure 21. Density of two shearing layers simulations with M = 1, at 
tkh . both without artificial conductivity. Left panel: Cubic spline kernel. 
Right panel: LIQ kernel. As previously shown the particle grouping of the 
CS kernel is abscent when using the LIQ kernel. Using the LIQ kernel also 
prevents the formation of a wide "gap" between the layers. 
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Figure 22. Density plots of shearing layers simulations with LIQ3 setup 
(see tables 3 and 4, M = 0.6), using different signal velocities. Left col- 
umn: t = Tkh- Right column: t = 2. Top two rows (label vi): use the 
signal velocity from eq. (5). Bottom two rows (label V2): use the signal 
velocity from eq. (6). lowres: 184 100 particles, highres: 552 100 particles. 



ii) might very well be a consequence of point i). We find that it is 
easy to actually lose substructure that was initially resolved in SPH 
because of superfluous energy diffusion (see Figs. 11, 20 and 22). 
When using SPH to solve very sensitive hydrodynamical problems 
(like the current shearing layers problem) one should therefore take 
good care when including the artificial conductivity and in selecting 
an appropriate signal velocity. We presented two new signal veloci- 
ties to that effect: v^^^^i (eq. (5)) and Vsig^^2 (^Q- (6))- Shearing layers 
results indicate that both these signal velocities lead to less energy 
diffusion in simulations. v^ig,i emerges as the best choice, giving 
better results in terms of KH rolls and requiring no significant extra 
computations. 

The combined effect of both the shock waves and the parti- 
cle disorder becomes less and less important as the time-scale of 
the problem itself (in this case: tkh) decreases. At he resolution 
of current galaxy formation simulations mixing is probably not im- 
portant. However, mixing could become crucial for next-generation 
simulations. Also, for the standard astrophysical application the ac- 
tual choice of AC sigal velocity will be less of a concern, if a con- 
cern at all, as is demonstrated by the host of successful problems 
tackled by the SPH codes of Rosswog & Price (2007) and Kawata 
et al. (2009). Note that a suitable signal velocity for simulations 
including gravity has not as of yet been presented. 
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Calculating is then straightforward. We give the expressions here 
for completeness. In two dimensions: 




(A9) 



Wr{u)du = 27r(|4 + f 4 + ^4 (AlO) 

+ f ^S + fx3')[ . (All) 

In three dimensions: 

Wr{u)du = 47T ^^Fxl - ^xi^ (A12) 

J' Wr(u)du = 4Tr(^jxl+jxt+jxl (A13) 

+ ^xt + ^x!)\' . (A14) 



APPENDIX A: LIQ KERNEL COEFFICIENTS 

The expressions for the parameters of the LIQ kernel (eq. (9)), ob- 
tained by solving the set of equations (lOa-lOf), are: 
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with Xs the LIQ kernel connection point. From this it is straightfor- 
ward to calculate the norm N by integrating over the volume: 
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